N = [2:65,70:5:200]';
x = load('tableData.dat');

Hnumerical = x(:,2) ./ 2;

Hintegral = N./8 .* (-2.*sqrt(N) + 4.*N  - 2 + 4.*log(2) - 2.*log(N) - 4.*log(2).*N - 2.*N.*log(sqrt(N)+1) + N.*log(N) + 2.*log(sqrt(N)+1));

H1 = N.^2 /2 - 1/2 * N.^(3/2) + 1/12 * N.^(1/2) - 1/240 * N.^(-1/2) + 1/10080 * N.^(-3/2) - 1/725760 * N.^(-5/2);

H2 = N.^2 / 8 * (1 - 2*log(2)) - 1/8 * N.*log(N) - 1/8*N*(1-2*log(2)) + 1/24 * log(N) - 1/12 * log(2) + N.^(-1) * (1/144 + 1/90 * log(2)) - 1/180 * N.^(-1).*log(N) - N.^(-2)*(7/6480 - 1/1260 * log(2)) + 1/2520 * N.^(-2).*log(N);

H3 = -N.^2/8 *(1+2*log(2)) + log(2)/4 * N + 1/6 * N.^(1/2) - (1/16+1/12* log(2)) - 1/20 * N.^(-1/2) + N.^(-1) * (1/48 + 1/90 * log(2)) + 23/5040 * N.^(-3/2) - N.^(-2)*(1/480 + 1/1260*log(2))-31/90720 *N.^(-5/2);


plot(N,Hnumerical,N,Hintegral - N.*log(N)/4);

%plot(N,Hnumerical,N,Hintegral,N,H1+H2+H3);

%plot(N,(Hintegral - Hnumerical)./Hnumerical); xlim([20 200]);
